# -*- coding: utf-8 -*-
"""
Created on 03 08 下午4:55 2022 
@Author : HHQUAN
@Email : 1075960398@qq.com
3516数据读取方法，自相关的方法有点问题
"""

import numpy as np
import matplotlib.pyplot as plt

# fid = open('/media/cl/cldata/3516data/2022_06_28/3516/4.3.4/rec.pcm', 'rb')
# tmp = fid.read()
# fid.close()
# data1 = np.frombuffer(tmp, dtype=np.int16, count=-1)
# sig = np.reshape(data1, newshape=(-1, 6))
# fid = open('/media/cl/cldata/3516data/2022_06_28/3516/4.3.4/out.pcm', 'rb')
# tmp = fid.read()
# fid.close()
# out = np.frombuffer(tmp, dtype=np.int16, count=-1)

fid = open('/home/cl/Downloads/2020下半年/speex相关/20220518测试/speexdsp-1.2rc3/release/build/data8_2/mic0.pcm', 'rb')
tmp = fid.read()
fid.close()
mic = np.frombuffer(tmp, dtype=np.int16, count=-1)

fid = open('/home/cl/Downloads/2020下半年/speex相关/20220518测试/speexdsp-1.2rc3/release/build/data8_2/out.pcm', 'rb')
tmp = fid.read()
fid.close()
out = np.frombuffer(tmp, dtype=np.int16, count=-1)

# plt.figure()
# plt.subplot(211)
# plt.plot(mic)
# plt.subplot(212)
# plt.plot(out)
# plt.show()

fs = 48000
fftLen = 960*2  # 长度对延时估计的影响
overlap = fftLen//2
frames = mic.size//overlap-1
win = np.hanning(fftLen).reshape((fftLen, ))
df = fs//fftLen

tdoa = np.zeros((frames, ), dtype=np.int)
eng = np.zeros((frames, ), dtype=np.float)
corrxy = np.zeros((fftLen, frames), dtype=np.float32)

for i in range(frames):
    block1 = np.fft.rfft(mic[i*overlap: i*overlap+fftLen]*win + 0.1*np.random.randn(fftLen))
    block2 = np.fft.rfft(out[i*overlap: i*overlap+fftLen]*win + 0.1*np.random.randn(fftLen))
    cross_spectrum = block1*np.conjugate(block2)/np.abs(block1*np.conjugate(block2))
    # fl, fh = 500//df, 8000//df
    # cross_spectrum2 = np.zeros_like(cross_spectrum)
    # for j in range(fl, fh, 1):
    #     cross_spectrum2[j] = cross_spectrum[j]
    xy = np.fft.ifftshift(np.fft.irfft(cross_spectrum))
    tdoa[i] = np.argmax(xy) - fftLen//2
    eng[i] = xy[np.argmax(xy)]
    corrxy[:, i] = xy # 这个是用来分析门限的参数

coxy2 = corrxy[fftLen//2-100:fftLen//2+100, :]
location = np.argmax(coxy2, axis=0)
plt.pcolormesh(coxy2)
plt.colorbar()
for n in range(location.size):
    if coxy2[location[n], n]>0.2: # 选择合理的门限可以有效减小野值，尽量不使用25毫米的间距
        plt.plot(n, location[n], '.k')
plt.show()




# L = 4800*10
# frames = sig.shape[0]//L
# maxn = np.zeros((frames, ))

# for i in range(frames):
#     x = sig[i*L:i*L+L, 0]
#     y = out[i*L:i*L+L]
#     xf1 = np.fft.rfft(x)
#     xf2 = np.fft.rfft(y)
#     cross = xf1*np.conj(xf2)/np.abs(xf1*np.conj(xf2))
#     xx1 = np.fft.ifftshift(np.fft.irfft(cross))
#     #
#     nn = np.argmax(xx1)
#     if xx1[nn]>0.3:
#         maxn[i] = nn
#     else:
#         maxn[i] = -1
#
# plt.figure()
# plt.plot(maxn)
# plt.plot(np.ones(maxn.shape)*L/2, ':r')
# plt.show()

